Matlab
code 2.6: Matlab file “Figure 2-11.m”
%--------------------------------------------------------------------
% This code can be used to generate Figure 2.11
%--------------------------------------------------------------------
clear all;
close all;
%% the theoretical curve of precession target
micro-Doppler characteristic when radar platform is vibrating
c = 3e8;
j = sqrt(-1);
fc = 10e9; % carrier frequency of transmitted
signal
v = 0; % translational velocity of target
cord = 1000*[300 100 500]; % coordinates of
local coordinate system's origin in the radar coordinate system
colo = [0 0 0]; % coordinates of radar in the
radar coordinate system
R0 = cord-colo;
n = R0/sqrt(sum(R0.^2)); % unit vector of
Line-of-Sight(LOS)
p1 = [0 0 0.5]; % scatterer on the apex of cone
p2 = [0.5 0 -0.5]; % scatterer on the bottom of
cone
tgt = [p1;p2];
alpha = atan(cord(2)/cord(1)); % azimuth angle
beita =
atan(cord(3)/(sqrt(cord(1)^2+cord(2)^2))); % elevation angle
alpha0 = pi/10; % azimuth angle of radar
platform vibrating axis
beita0 = pi/3; % elevation angle of radar
platform vibrating axis
dr = 0.01; % vibration amplitude of radar platform
fr = 8; % vibration frequency of radar platform
ae = pi*[1/3 1/4 1/5]; % Euler angle
ws = [19.2382 -11.1072 22.2144]; % angular
velocity of spinning in the local coordinate system
omegas = sqrt(sum(ws.^2));
% omegas = 10*pi;
Ts = 0.2; % spinning period
wc = [7.6756 -2.3930 9.6577]; % angular
velocity of coning in the local coordinate system
omegac = sqrt(sum(wc.^2));
% omegac = 4*pi;
Tc = 0.5; % coning period
t = 2;% radar illumimated time
ri1 = [cos(ae(3)) sin(ae(3)) 0;-sin(ae(3))
cos(ae(3)) 0;0 0 1];
ri2 = [1 0 0;0 cos(ae(2)) sin(ae(2));0
-sin(ae(2)) cos(ae(2))];
ri3 = [cos(ae(1)) sin(ae(1)) 0;-sin(ae(1))
cos(ae(1)) 0;0 0 1];
ri = ri1*ri2*ri3; % initial rotation matrix
ws = ri*ws'; % angular velocity of spinning in
the reference coordinate system
% omega = sqrt(sum(w.^2));
wse = ws/omegas; % unit vector of spinning
wsr = [0 -wse(3) wse(2);wse(3) 0
-wse(1);-wse(2) wse(1) 0]; % skew symmetric matrix of spinning
wc = ri*wc'; % angular velocity of coning in
the reference coordinate system
wce = wc/omegac; % unit vector of coning
wcr = [0 -wce(3) wce(2);wce(3) 0
-wce(1);-wce(2) wce(1) 0]; % skew symmetric matrix of coning
r0a = tgt(1,:)-colo;
r0b = tgt(2,:)-colo;
r0ar = ri*r0a'; % spinning scatterer in the
reference coordinate system
r0br = ri*r0b'; % coning scatterer in the
reference coordinate system
a1 =
(2*fc/c)*r0br'*(omegac*wcr^2+omegac*wcr^2*wsr^2)'*n';
a2 =
(2*fc/c)*r0br'*(omegac*wcr+omegac*wcr*wsr^2)'*n';
a3 =
(2*fc/c)*r0br'*(omegac*wcr^2*wsr+omegas*wcr*wsr^2)'*n';
a4 = (2*fc/c)*r0br'*(omegas*wcr*wsr-omegac*wcr^2*wsr^2)'*n';
a5 =
(2*fc/c)*r0br'*(omegac*wcr*wsr-omegas*wcr^2*wsr^2)'*n';
a6 =
(2*fc/c)*r0br'*(-omegac*wcr*wsr^2-omegas*wcr^2*wsr)'*n';
a7 =
(2*fc/c)*r0br'*(omegas*wsr^2+omegas*wcr^2*wsr^2)'*n';
a8 = (2*fc/c)*r0br'*(omegas*wsr+omegas*wcr^2*wsr)'*n';
prf = 4000; % pulse repetition frequency
pri = 1/prf; % pulse repetition interval
dt = 0:pri:t-pri; % time sampling interval
m = length(dt);
fmd1 = zeros(2,m); % theoretical micro-Doppler
frequency
fmd2 = zeros(2,m);
for i = 1:m
fmd1(1,i) =
a1*sin(omegac*dt(i))+a2*cos(omegac*dt(i))+a3*sin(omegas*dt(i))*sin(omegac*dt(i))...
+a4*sin(omegac*dt(i))*cos(omegas*dt(i))+a5*cos(omegac*dt(i))*sin(omegas*dt(i))...
+a6*cos(omegac*dt(i))*cos(omegas*dt(i))+a7*sin(omegas*dt(i))+a8*cos(omegas*dt(i))...
-(4*pi*fc*fr*dr/c)*(cos(alpha-alpha0)*cos(beita)*cos(beita0)+sin(beita)*sin(beita0))*cos(2*pi*fr*dt(i));
fmd1(2,i) =
(2*fc*omegac/c).*((wcr^2.*sin(omegac*dt(i))+wcr.*cos(omegac*dt(i)))*ri*r0ar)'*n'...
-(4*pi*fc*fr*dr/c)*(cos(alpha-alpha0)*cos(beita)*cos(beita0)+sin(beita)*sin(beita0))*cos(2*pi*fr*dt(i));
fmd2(1,i) =
a1*sin(omegac*dt(i))+a2*cos(omegac*dt(i))+a3*sin(omegas*dt(i))*sin(omegac*dt(i))...
+a4*sin(omegac*dt(i))*cos(omegas*dt(i))+a5*cos(omegac*dt(i))*sin(omegas*dt(i))...
+a6*cos(omegac*dt(i))*cos(omegas*dt(i))+a7*sin(omegas*dt(i))+a8*cos(omegas*dt(i));
fmd2(2,i) =
(2*fc*omegac/c).*((wcr^2.*sin(omegac*dt(i))+wcr.*cos(omegac*dt(i)))*ri*r0ar)'*n';
end
figure(1)
plot(dt,fmd1,'r') % when radar platform is
vibrating
hold on
plot(dt,fmd2,'b') % when radar platform is
stationary
xlabel('Time (s)')
ylabel('Frequency (Hz)')
axis([0,2,-1200,1200])
%% the time-frequency analysis result of
precession target micro-Doppler characteristic when radar platform is vibrating
r = zeros(length(tgt(:,1)),length(dt)); %
distance between the scatterers and radar
for i = 1:m
for n
= 1:length(tgt(:,1))
if n == 1
r(n,i) = sqrt((R0'+v*t+((eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))*r0ar)...
-(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)]))'...
*(R0'+v*t+((eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))*r0ar)...
-(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)])));
else
r(n,i) =
sqrt((R0'+v*t+(eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))...
*(eye(3)+wsr*sin(omegas*dt(i))+wsr^2*(1-cos(omegas*dt(i))))*r0br...
-(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)]))'...
*(R0'+v*t+(eye(3)+wcr*sin(omegac*dt(i))+wcr^2*(1-cos(omegac*dt(i))))...
*(eye(3)+wsr*sin(omegas*dt(i))+wsr^2*(1-cos(omegas*dt(i))))*r0br...
-(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)])));
end
end
end
s = zeros(length(dt),length(tgt(:,1))); % echo
signal matrix
for i = 1:length(tgt(:,1))
s(:,i) = exp(j*2*pi*fc*2*r(i,:)'/c);
end
st = sum(s,2);
N = 200; % number of Gabor coefficients in time
Q = 50; % degree of oversampling
tfr = tfrgabor(st.',N,Q); % Gabor transform
tfr_r = fftshift(tfr,1);
figure(2)
contour(linspace(min(dt),max(dt),length(tfr_r(1,:))),linspace(-1000,1000,length(tfr_r(:,1))),tfr_r,70)
xlabel('Time (s)')
ylabel('Frequency (Hz)')
axis([0,2,-800,800])